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ABSTRACT 

Coincident detections of electromagnetic (EM) and gravitational wave (GW) signatures from coalescence 
events of supermassive black holes are the next observational grand challenge. Such detections will provide 
the means to study cosmological evolution and accretion processes associated with these gargantuan compact 
objects. More generally, the observations will enable testing general relativity in the strong, nonlinear regime 
and will provide independent cosmological measurements to high precision. Understanding the conditions 
under which coincidences of EM and GW signatures arise during supermassive black hole mergers is therefore 
of paramount importance. As an essential step towards this goal, we present results from the first fully general 
relativistic, hydrodynamical study of the late inspiral and merger of equal-mass, spinning supermassive black 
hole binaries in a gas cloud. We find that variable EM signatures correlated with GWs can arise in merging 
systems as a consequence of shocks and accretion combined with the effect of relativistic beaming. The most 
striking EM variability is observed for systems where spins are aligned with the orbital axis and where orbiting 
black holes form a stable set of density wakes, but all systems exhibit some characteristic signatures that can 
be utilized in searches for EM counterparts. In the case of the most massive binaries observable by the Laser 
Interferometer Space Antenna, calculated luminosities imply that they may be identified by EM searches to 
jsil, while lower mass systems and binaries immersed in low density ambient gas can only be detected in the 
local universe. 

Subject headings: accretion - black hole physics - gravitational waves - hydrodynamics - relativity - shock 
waves - X-rays: bursts 
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1. INTRODUCTION 

Galactic mergers are expected to be the major route 
for the formation of binary black holes (BBHs), in 
agreement with predict ions of hierarchical mo d els of 
structure formation ([Haehnelt and Kauffmannl 120021: 
IVolonteri. Haardt. and Madaul 12003 ah and with obser- 
vations showing that the majority of galaxies harbor 
supermassive black hole s (SMB Hs ') in their cente rs (e.g. , 
Kormendv and Richstonel J 1 9951: iRichstone et ail 119981: 



Peterson and WandeW2000t Ferrarese and Fordll2005l) . some 



fraction of newly formed galaxies are thus expected to host 
binary systems of SMBHs, whose inspiral and merger lie in 
the "sweet spot" of the future Laser Interferometer Space 
Antenna (LISA). These cataclysmic cosmic events produce 
copious amounts of gravitational radiation (~ 10 57 ergs _1 ) 
during the final stages of merger, radiation which is imprinted 
with detailed information about the inspiral and merger 
of the binary as well as the state (mass, spin and kick) 
of the resulting black hole (BH). Detecting and charac- 
terizing these binary mergers will thus be an important 
resource for understanding formatio n mechanisms and the 
cosmological evolution of SMBHs ( Volonteri et al.l l2003bl 
I2004t iMicic et al.l [20071: ISesana et al.1 120071) . A synergistic 
detection of electromagnetic (EM) and gravitational wave 
(GW) radiation from merging SMBHs h as been defined as 
the next observational grand challenge dKocsis et al.l 120061, 
for example). A combination of EM and GW signatures 
will provide independent measurements of redshift and 
luminosity distance, allowing cosmological measurements 
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at a precision of ~ 10%, limited mainly by uncertaintie s 
due to weak gravitational lensing dHughes and Holzll2003l) . 
Similarly, EM+GW observations can be used to test whether 
gravitons travel at the speed of light, as required by general 
relativity, and in such a way test one of the fundamental 
principles of the theory. Furthermore, these multi-messenger 
observations will provide new insight on accretion processes 
associated with SMBH binary systems in the final stages of 
their evolution. 

One of the outstanding astrophysical questions with di- 
rect bearing on the feasibility of EM+GW detections is re- 
garding the physical properties of the gaseous environment 
surrounding a binary before and during coalescence. Ob- 
servationally, a significant sample of these objects is yet to 
be attained, and finding them in EM searches is a chal- 
lenging task (iBogdanovic et al.ll2009h . Consequently, most 
of the information about these systems has to be drawn 
from theoretical perspectives. Non-relativistic hydrodynam- 
ical simulations have, over the past several years, signifi- 
cantly contributed to our understanding of the evolution of 
BH pairs an d interstellar gas, both during and after the galac- 
tic mergers (Kazantzidis et al. 2005; Armitage and Natarajan 



2002; 



2007 



Escala et alj|2004 l2005t iDotti et al.ll2007t iM aver et al 
Colpi et alj|2007t iMacFadven and Milosavlievicll2008l 
Havasaki et al.l 120081: ICuadra et al.l 120091) . However, simu 
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lations spanning the entire dynamical range, from galactic 
merger scales (~ 10 2 kpc) to binary coalescences (< 10~ 2 pc), 
are still prohibitively computationally expensive. As a con- 
sequence, non-relativistic simulations stop at binary separa- 
tions of order 1 pc while fully general relativistic simulations 
are possible only at separations of order 10~ 5 pc. Thus the 
properties and structure of accretion flows around binaries 
are uncertain. For instance, the presence of gas on larger 
scales in the aftermath of a gas rich galactic merger may 
not guarantee an abundant supply of gas for accretion once 
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a gravitationally bound binary is formed. This is because 
the binary torques can evacuate most of the surrounding gas, 
thus preven ting any significant accretion on e ither member of 
the binary dMilosavlievic and Phinnevll2005l) . The scenario 
in which binary torques clear a central low density region 
is commonly described in the literature as the circumbinary 
disk. Alternatively, it is also plausible that if the surround- 
ing gas is sufficiently hot and tenuous, the binary may find 
itself engulfed in a radiatively inefficient, turbulent flow all 
the way through coalescence. Such conditions could arise 
in gas deficient mergers and are indeed expected to exist in 
nuclear regions of some low luminosity active galactic nuclei 
(AGNs:l5uataertll999HPtak et al.l2004HNemmen et d]2006t 
lElitzur and Holl2009l for example). We refer to this scenario 
as the gas cloud and note that if binaries indeed do exist in 
radiatively inefficient flows, then the circumbinary disk and 
gas cloud scenarios effectively bracket the range of physical 
situations in which pre-coalescence binaries may be found in 
centers of galaxies. 

In addition to EM signatures and accretion, a gaseous envi- 
ronment could potentially have a profound effect on the dy- 
namics of the BHs. For instance, accretion torques in gas 
rich mergers could force the BHs to have a "preferential" 
spin orientation, aligning the spins of the B Hs with the angu- 
lar m omentum of the large-scale gas disk dBogdanovic et alj 
l2007h . Since numerical relativity (NR) simulations have 
established that the magnitude of the gravitational recoil 
on the final B H depends on the sp i n orientation of the 
merging BHs dHerrmann et al.l l2007at iKoppitz et all 120071: 



Campanelli et al. 2007a b; Gonzal ezet al.l 120071: iBaker et al] 
20071: ISchnittman and Buonannoll2007l) . the presence or ab- 
sence of gas in the vicinity of the BBH can have direct im- 
plications for the magnitude of the kick inflicted on the final 
BH. 

In this work, we build upon the framework of larger- 
scale simulations as well as the initial relativistic calcu- 
latio ns that investigated b oth the dynamics of test parti- 
cles (Ivan Meter et all 120091) and the evolution of EM fields 
dPalenzuela et alJl2009allbF in the gravitational potential of a 
coalescing binary. We use fully general relativistic numer- 
ical hydrodynamics simulations to follow the interaction of 
a BH binary in a gaseous environment through coalescence. 
We focus on the final stages of the binary evolution (last few 
orbits and merger) and consider only equal-mass SMBH bi- 
naries surrounded by a hot and tenuous gas cloud. The main 
objective of this work is to characterize the EM and GW sig- 
natures that arise during coalescence 4 . 

The remainder of this paper is organized as follows: The 
computational methodology used is described in ^2j followed 
by the initial conditions used in the simulations in |3] We 
describe the general gas dynamics in $4] In f5] we present 
a discussion of the EM and GW signatures, followed by our 
conclusions in $6] 

2. COMPUTATIONAL METHODOLOGY 

The results in this paper were obtained with the new ver- 
sion of the MayaKranc code of the NR group at Georgia 
Tech. The new code is an enhanced version of the code that 
was used primarily for studi es involving vacuum spacetimes 
conta i ning BH singularities dWashik et a l. 2008; Hinder et al 
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2009). As with the previous code, the new MayaKranc 
code is base d on the BSSN formulation and the moving punc- 
ture method dCampanelli et al.ll2006aUBaker et all2006l). The 
code source is generated by the package Kranc dHusa et all 
120061) . which produces a set of thorns t hat work under the 
Cactus in frastructure dAllen et al.ll 19991) and Carpet mesh 
refinement (ISchnetter et alj|2004l) . The new main feature in 
the MayaKranc code is the inclusion of general relativistic 
hydrodynamics. The hydrodynamics code in MayaKranc 
is a modified version of the public version of the WHISKY 
hydrodynamics code developed by the Europe an Union Net- 
work on Sources of Gravitational Radiation ( Whis ky^ webl : 
iBaiotti et al.ll2003l) . Unless otherwise specified, we use ge- 
ometrized units where G = c = 1 and normalized to the mass 
of the system M, using the metric signature (-,+,+,+). 

MayaKranc assumes a perfect fluid with energy momen- 
tum tensor 

TV = phu^ l u v + Pg til ' (1) 

with p the rest mass density, h= l+e + P/p the enthalpy, P 
the pressure, e the internal energy per unit mass, and u^ the 
4-velocity of the fluid. The 3-velocity of the fluid is given by 
v' = (u'/u° + f3')/a, with a and j3' the lapse function and shift 
vector, respectively. The quantities p, v' and e are considered 
primitive variables. The pressure P and the Lorentz factor 
W = au° are viewed as auxiliary variables. The pressure P 
for all the simulations in the present work is calculated from 
the T-law equation of state P = pe(T-l) with T = 5/3. The 
Lorentz factor is obtained from u^u^ = — 1. 



WHISKY uses the Valencia formulati on dMartfet al.ll 199 It 
iBanvuls et al.lll997t llbanez et al.l 1200 ll) of numerical hydro- 
dynamics. That is, the evolution equations are cast into a set 
of conservation equations of the form 



d F°(w) + d i F i (w) = S(w), 



(2) 



where w = (p, v ! , e) is the vector of primitive variables and F° = 
(D,S',t) the corresponding vector of conservative variables 
defined as: 



D = ^pW 

s^^pimy 

T=y^(phW 2 -P)-D 



(3a) 
(3b) 
(3c) 



4 Animations of the simulations discussed in this paper can be found at 

http : // www . era . gatech . edu/numrel / paper s /BBH_Gas Cloud 



with 7 the determinant of the spatial metric. The vectors F' 
and S are the fluxes and source terms, respecti vely. For details 
see lBanvuls et alj d 1997b : llbanez et al.l dioOll) . 

Our code, MayaKranc, modifies the Whisky code in 
several important ways. In order to improve efficiency as well 
as simplify the interface between the hydrodynamics and ge- 
ometry sectors, we implemented a new construction of the 
stress-energy tensor and the corresponding sources ap- 
pearing in the BSSN equations. Furthermore, in the regions 
wi thin the apparent ho rizons (AHs), we impose, as suggested 
by iFaber et alj d2007l) . the dust limit of the hydrodynamics 
equations when unphysical data appears during the evolution. 
We also modified the atmosphere treatment used to model 
vacuum regions. A filter is in place in the atmosphere domains 
to avoid spurious increases of fluid densities above the atmo- 
spheric threshold limit during temporal updates. The overall 
structure of the Whisky code has also slightly changed to 
be able to work seamlessly with the rest of the MayaKranc 
code. 

For BBH initial data, the MayaKranc c ode uses the 
shtiiffUNCTURES spectral code developed by lAnsorg et alj 
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d2004l) . We modified this solver to include matter fields. Since 
the present work focuses on matter fields initially at rest, 
the modifications to the 2PUNCTURES code involved only 
the Hamiltonian constraint. The solutions to the momentum 
constraint remained unchanged; specifically the Bowen-York 
type solutions continue to be applicable dBowen and Yorkl 
198(f). As with the vacuum case, we construct initial data 
for a BBH in a gaseous environment assuming both confor- 
mal flatness and a vanishing trace of the extrinsic curvature. 
The Hamiltonian constraint thus reads 

V 2 V>+ l^AijA'i = -2n^p = -27rV>" 3 po , (4) 
8 

where p Q = n^n^T^ = phW 2 -P with the unit normal to 
the constant time hypersurfaces. In Eq. [4] conformal quanti- 
ties have tildes, Ay is the traceless part of the extrinsic curva- 
ture, and ip is the conformal factor such that the spatial metric 
transforms like 7y = ip 4 T]ij with rjtj the flat metric. The con- 
formal rescaling p a = tp s p is necessary in order to guar antee 
the ex istence of a solution to the constraint equation dYorkl 
119791) . It is important to emphasize that the freely specifiable 
matter field is p Q . We construct p from p = ip„ ac po in which 
ipvac is the solution to Eq. [4] in the absence of matter. Given 
p constructed in this fashion, the solution to the Hamiltonian 
constraint yields a different ip, and thus the new p Q does not 
exactly correspond to the primitive matter fields, p , used in 
setting up the conformal factor. 

The MayaKranc code also includes infrastructure to an- 
alyze the data produced by the simulations. This infrastruc- 
ture consists of a set of analysis tools or modules to construct 
gravitational waveforms, estimating BH kicks and spins, and 
tracking AHs. With the inclusion of hydrodynamics, more 
modules were needed. Of particular relevance to the present 
work, we developed tools to monitor conservation of total rest 
mass in the system, to compute fluxes of matter both through 
AHs (accretion) and in the form of outflows from the imme- 
diate vicinity of the BHs, and to calculate the integrated lumi- 
nosity of the emitted light. 

We tested the implementation of hydrodynamics in 
MayaKranc through a suite of standard tests. One of the 
important tests includes verifying the conservation of bary- 
onic mass. While baryonic mass should be conserved ev- 
erywhere in the spacetime, we consider here the conservation 
within sub-volumes of the computational domain, which also 
aids us in testing in which sub-volu mes conservation ma y be 
problematic. That is, integrating (see lBanvuls et all 1997b 

d,D + d i [D(av i -p i )]=0 (5) 
over a constant volume V with boundary dV yields 

M v + [ D(av i -(3 i )r i dS = 0, (6) 

JdV 

with f the outgoing unit vector to V and 

M v = f Dd 3 x. (7) 
Jv 

Time integration of Eq.|6]yields 

M(t) = Mv(t) + M dv (t) = constant. (8) 

where 

M 9V (t)= ( dt' ( D(a^-0)ndS. (9) 
Jo Js h 



On a given step, we compute Ai from Eq.[8]to check how well 
it remains constant. 

As a test of the code, we present here th e evolution of 
a Tolman-Oppenheimer-Volkoff (TOV) star (iTolmanl 119391: 
lOppenheimer and VolkofdU"939t iMisner eTa!lll97l " of mass 
I AM, where M is the arbitrary mass scale in the simulation. 
A TOV star is a static, spherically symmetric solution to the 
Einstein equations for a spacetime with a perfect fluid. As 
such, the stability of an evolved TOV is a standard test in 
numerical relativity with hydrodynamics. The initial data is 
created by solving the TOV equations, assuming a polytropic 
equation of state P= np r with k = 100 and T = 5/3. The com- 
putational grid for these test runs consisted of 5 refinement 
levels extending to an outer boundary of 512M, the finest of 
which lay within the star. We evolved this system for three 
resolutions where we chose the resolutions on the finest re- 
finement level as M/8, M/11.25, and M/16. In the top panel 
of Fig. [T] we plot the errors in baryonic mass conservation, 
SM/Mq = (M-Mo)/Mo, as a function of time with M 
computed from Eq. [8] and A4o = Ai(t = 0). In the bottom 
panel of Fig. Q] we show the corresponding errors in the cen- 
tral density, p mdx , over 500M of evolution time. The central 
density stays within an acceptable 1.5%. The baryonic mass 
conservation on the other hand shows a linear increase of er- 
rors for the coarsest level of resolution of M/8. This error 
growth arises from the refinement boundaries in the vacuum 
regions, which are modeled as a low density atmosphere. As 
each boundary comes into causal contact with the surface of 
the star, the density rises above atmospheric levels despite our 
atmospheric filters, adding a mass violation. While this ef- 
fect is significant in the coarsest run, the errors remain below 
0.2% in the other, higher resolution runs. Improvements in 
atmosphere handling are thus needed in our code for systems 
involving vacuum regions (e.g. BH + neutron star mergers). 
Later on, in Section |4] we show that errors in mass conser- 
vation do not affect the results presented in this paper since 
the relevant dynamics remain completely contained within the 
gas cloud. 
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FIG. 1. — Top panel shows the errors in baryonic mass conservation as a 
function of time from the evolution of a TOV at different numerical resolu- 
tions. Bottom panel shows the corresponding errors in the central density. 

3. INITIAL DATA 

With the uncertainties in our understanding of accretion 
flow properties around binaries, it is not possible to uniquely 
define initial conditions for simulations of BBHs and gas on 
very small scales (~ 10~ 5 pc), such as those considered in the 
present work. Nevertheless, the physical properties of the gas 
in the vicinity of the BHs are likely to be bracketed by two 
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characteristic scenarios depending on the balance of heating 
and cooling mechanisms in the accretion flow. In the first 
scenario, a hot and turbulent accretion flow forms as a con- 
sequence of relatively inefficient cooling processes; the BBH 
is then immersed in a pressure-supported, geometrically thick 
torus with a scale height H comparable to its radius R. Here- 
after, we refer to this scenario as the gas cloud model. In 
the second scenario, relatively efficient cooling processes re- 
sult in the gas settling into a rotationally-supported, geometri- 
cally thin (H < R) accretion disk around the binary. This sce- 
nario is commonly referred to as the circumbinary disk model 
and has been previously stud ied by means of non-rela t ivistic 
hydrodynamic simulations (lArmitage and Nataraianl 
Milosavljevic and Phinnev 2005; Havasaki et al. 2007, 



2002; 



2008 



MacFadven and Milosavlievidl2008t Icuadra et alj|2009t) . In 



this paper, we consider the gas cloud scenario and defer the 
numerical investigation of the circumbinary disk with a fully 
relativistic treatment to a future study. 

Though the majority of BBHs in the universe are expected 
to be unequal-mass binaries, we will consider for computa- 
tional accessibility and as a first step only equal-mass bina- 
ries where we can exploit symmetries of the system. Our 
simulations thus consist of equal-mass binaries with total 
mass M and initial separation of R = 8M « 10~ 5 M7pc, where 
M-j =M/10 7 M Q . This choice of mass scaling yields BBHs 
which can be detected in the LISA band during the plunge 
and coalescence and therefore in a regime where modeling 
is only accessible via NR. As these systems are among the 
higher mass BBH systems in LISA'S sensitivity window, they 
are relatively luminous in gravitational radiation and, if asso- 
ciated with an EM signature, likely to be among the most EM 
luminous binaries that LISA can detect. 

The initial data for our simulations consist of a BBH im- 
mersed in a gas cloud with a radius of 6QM. The radius of the 
cloud is selected arbitrarily and in such a way that it entirely 
encloses the binary orbit. The gas cloud is initially static with 
a Gaussian density profile 



p= p Q e 2„i 



(10) 



where a = 10.83M and p is the central density in the cloud. 
The a was chosen to avoid non-atmosphere gas densities from 
reaching the outer boundary initially, though subsequent sim- 
ulations show larger a do not change the results qualitatively. 
The value of p Q adopted in our simulations is obtained by 
extrapolating the results of larger scale simulations that fol- 
low the evolution of BBHs in gaseous environment in the af- 
termath of galactic mergers. These results suggest that the 
amount of gas which remains strongly bound to the binary on 
scales of < 10~ 2 pc can be of the or der of 1% of the total mass 
M of the binary dColpi et al.ll2007l) . This implies a mean gas 
density of 



0.01M 
(0.01 pc) 3 



10" n M 7 gcm" 



(11) 



Note that the internal units in the MayaKranc code are 
given in terms of the total mass of the binary, M. For vac- 
uum simulations, this implies that the results from a BBH 
simulation with mass M can be scaled to arbitrary BH phys- 
ical masses. For non-vacuum simulations such as those in 
the present work, the BH mass also determines the scaling 
of hydrodynamical variables (e.g,. density, pressure, internal 
energy, etc.). In particular, density scales as M~ 2 . Thus, in 
non-vacuum simulations scaling with M is not arbitrary and 



the mass parameter should be chosen in such a way as to re- 
flect a plausible range of densities. We set the initial central 
density of the gas cloud to p a = 7 x 10~ 12 M7 2 gcm~ 3 , a value 
consistent with the mean density estimate in Eq. QT| The to- 
tal rest mass of the gas cloud in the computational domain is, 
initially, ~ 10~ 4 M7M Q , about 11 orders of magnitude lower 
than the BH masses. For computational facility, we surround 
the gas cloud with a uniform, low density atmosphere of den- 
sity 10~ 5 p o . 

We use a polytropic equation of state P = k p T to construct 
the initial data; the temperature of the gas is then given by 



m„c 2 P m„c 2 



tfi p 



(12) 



In this gas cloud model, the thermal energy of the gas is com- 
parable to its gravitational potential energy. Thus, the velocity 
of the fluid in the vicinity of the BHs is 



v'th : 



l3k B T 



:0.35c 




(13) 



where T is the temperature of the gas. It follows that the 
gas in the vicinity of the binary is a high temperature plasma, 
T ~ 10 12 K, with thermal velocities comparable to the binary 
orbital speed. We then estimate the adiabatic constant k in the 
polytropic equation of state by evaluating Eq.Q~2]at the center 
of the cloud, where T = 10 12 K and p = 7 x lO^gcm"^ 2 , 
and obtain re = 2.51 x 10 8 (cm 3 /g) r ~ 1 , where T = 5/3. For 
these values and an ideal fluid equation of state, we find the 
speed of sound at the center of the cloud, 



1 1 fdP P dP 
h\dp p 2 de 



(T-l)(e + P/p) 



0.28 c. 



(14) 

Because the speed of sound is comparable though smaller 
than the thermal velocity and the orbital velocity of the binary, 
it follows that shocks can be expected to develop in vicinity 
of the BHs. 

The gas cloud surrounding the BBH is not in hydrostatic 
equilibrium. Namely, if the cloud were placed in a poten- 
tial well of a single BH, the top layers of the cloud would 
be gravitationally unbound and the gas cloud would expand, 
doubling its volume over approximately 130M of evolution. 
However, a condition for hydrostatic equilibrium in the case 
of an isolated BH is of limited utility for the binary scenario 
considered here. In a realistic case, the binary torques will act 
to unbind some fraction of gas and, as a result, the outer layers 
of the cloud in the binary scenario can be expected to be hot- 
ter (i.e., less gravitationally bound) compared to the cloud in 
equilibrium around a single BH. The simple initial conditions 
we choose emulate this behavior. Note however that in or- 
der to achieve more realistic initial conditions one should, in 
principle, evolve the cloud in the potential of a rotating BBH 
long before the inspiral and merger until quasi-steady state is 
achieved. Such relaxation is nontrivial and computationally 
expensive in the case of a general relativistic system. Instead, 
we start the evolution of the gas cloud in the initially frozen, 
non-rotating gravitational potential of the BBH for 32M of 
evolution time, or four times the crossing-time of the system. 
During this phase, the cloud passively evolves to form a gravi- 
tationally bound core with two density peaks, each associated 
with one of the BHs. After this period, we "unfreeze" the bi- 
nary and evolve it on its orbit for an additional ~ 60M (half 
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orbit), during which the gas adjusts to the dynamics of the 
binary. Any transients that develop as a consequence of this 
initial relaxation are the artifact of imperfect initial conditions 
and are not considered as a true signature of the system. 

The initial distribution of velocity vectors of fluid elements 
is isotropic about the center of mass of the system. Given the 
high thermal velocities, it follows that in such a flow, the ra- 
dial inflow speed of the gas will always be sufficiently high 
that the binary torques will be incapable of evacuating all 
of the gas and creating a low density "hole" in the center of 
the cloud, a property that is a distinguishing feature from the 
circumbinary disk model. The only mechanisms for heating 
and cooling of gas in our simulations are adiabatic compres- 
sion and expansion in the potential of the binary. We do not 
account for heating of the gas by the active galactic nuclei 
(AGN) nor for cooling by emission of radiation. 

The initial orbital parameters for the BBHs are calculated 
as though the binary were in vacuum. This is justified since 
the total mass in the gas cloud is negligibly small compared 
to the mass of the BHs. We calculate the initial momentum 
required for quasi-circular orbits from post-Newtonian (PN) 
equations. For the non-sp inning BHs, we evolve the 3rd order 
PN equations of motion dBuonanno et al.ll2006l) from a sep- 
aration of 40M to the desired separation. We also consider 
binaries with spinning BHs where spins are parallel and ei- 
ther aligned or anti-aligned with the orbital angular momen- 
tum. In this case we calculate the initial momentum of the 
BHs fro m the 3rd-ord er PN angular momentum found in Eq. 
(4.7) of (lKidderill995h . assuming purely tangential momenta. 
Table[T]lists parameters for the BBH systems discussed in this 
paper. For each BH, labeled 1 and 2,a\^/m is the dimension- 
less spin parameter, where m is the mass of each BH set in 
such a way that m/M = 1/2, with M the total mass of the bi- 
nary. As mentioned before, the BHs are placed at a coordinate 
distance 8M along the x-axis, with BHi located at x = AM and 
BH2 at x = -AM. In TableQ] P x and P y are the components of 
linear momentum for BHi ■ Finally, m, is the irreducible mass 
for the individual BHs computed from the area of their AH, 
and Madm is the total initial energy of the system. 

The simulations are run on a grid whose outer boundary is 
located at 3 17.44 M, with 9 levels of mesh refinement. The 

4 coarsest levels are fixed about the center of mass while the 

5 finest levels are constructed around and free to follow the 
BHs. Forruns GO, Gl, andG3, the finest resolution is M/51.6 
and the refinement boundaries are chosen such that the 5 finest 
have radii of 32 grid points while the 4 coarsest have radii of 
64 grid points. The mesh setup for G2 was altered somewhat 
as higher spins require a higher resolution at the puncture. The 
physical refinement boundaries of G2 remain stationary while 
the resolution is increased for a finest resolution of M/67.7. 
For all runs, we halve the evolved computational domain by 
using the reflection symmetry about the xy-plane. For runs 
GO, Gl, and G2, we halve the domain again using rotational 
symmetry. 



4. DYNAMICS OF THE BINARY AND GAS 

We now describe the general features and dynamics of the 
binaries and gas in our modeled systems. BBHs with an ini- 
tial separation of 8M evolve for approximately 3 orbits be- 
fore they plunge and finally merge. Since the mass of the gas 
cloud is small compared to the BBH, the binary dynamics are 
practically indistinguishable from the equivalent situation in 



Run 


ai/m 


a 2 /m 


P*/M 


py/M 


m,/M 


Madm/M 


GO 








-2.0902 x 10 -3 


0.11237 


0.5000 


0.9878 


Gl 


+0.4 


+0.4 





0.10862 


0.4893 


0.9875 


G2 


+0.6 


+0.6 





0.10677 


0.4736 


0.9874 


G3 


+0.4 


-0.4 





0.11237 


0.4893 


0.9878 



TABLE 1 
BBH INITIAL PARAMETERS. 



vacuum. The only differences among the cases we consid- 
ered a re due to the "hang-up" or delay induced by the BH 
spins (ICampanelli et alJl2006bl) . The excess of angular mo- 
mentum from the BH's spins result in runs Gl and G2 taking 
longer to merge compared to runs GO and G3. In run G2 for 
example, the hang-up extends the inspiral to about 5 orbits 
before merger. 

In Figures |2] and [3] we show snapshots of the baryonic rest 
mass density for the cases G2 and G3. The distinct features 
that arise in the gas during the binary's evolution are the den- 
sity waves that develop behind the inspiraling BHs and a bar 
of high density gas connecting the two BHs. This pattern of 
density distribution of gas is not unique to relativistic systems 
and was previously noticed in simulations of binary systems 
interacting with gas on a variety of scales, including binary 
galactic nuclei and TTauri stars. In order to examine whether 
the gas dynamics driven by the binary give rise to shocks, we 
calculate the Mach number, (3 S = v/c s , by comparing the ve- 
locity of the fluid element to the speed of sound measured 
locally in the simulation 5 . Figures|4]and|5]show sequences of 
snapshots with the distribution of the Mach number value for 
the gas in the plane of the binary of runs G2 and G3. In all 
of the panels, only transonic and supersonic parts of the flow 
were plotted, where (3 S > 1 . Early in the inspiral, the shocks 
are confined to the density wakes expanding externally to the 
binary orbit. Later in the inspiral, the bar of high density gas 
forms between the BHs, and its expansion gives rise to an- 
other shock region emanating from the central part of the bi- 
nary orbit. As the BBH evolves closer to merger, the central 
high density region decreases in size, being replaced after the 
merger by a steadier inflow in the immediate vicinity of the 
final BH. 

Here we briefly make a note of the degree to which bary- 
onic mass is conserved in our simulated binary systems, as 
discussed in the context of a TOV star at the end of ^2] We 
monitor the conservation of mass in two separate volumes: (1) 
the volume r < 8M excluding the regions inside the AHs and 
(2) the volume 8M < r < 20 M, external to the binary orbit. 
This allows us to calculate the total error in mass conserva- 
tion over the whole computational domain but also to more 
closely monitor the inner region containing the binary. From 
run GO we measure that the baryonic mass is conserved in 
both volumes at the level of < 1% over the duration of the 
simulation. Due to similarities in computational setup and the 
physical systems, similar mass conservation is expected in the 
other runs. 

5 . ELECTROMAGNETIC AND GRAVITATIONAL WAVE 
SIGNATURES 

In this section, we discuss the EM and GW signatures ex- 
pected to arise from the binary dynamics and accretion flows 

5 Though fi s > 1 is necessary, it is not a sufficient condition to locate 
shocks (e.g., cold fronts). Given the long cooling timescale of the hot, diffuse 
gas, however, discontinuities in the Mach number within our simulations do 
select regions where shocks are likely to form. 
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FIG. 2. — Snapshots of the baryonic rest mass density, p [gem 3 ], in the orbital plane of the binary shown for ran G2. From upper left to lower right the 
snapshots correspond to t = 7.6 X 10 2 Mi s, 1 .3 X 10 3 M7 s, 1 .9 X 10 3 M7 s, 2.5 X 10 3 M7 s, 3.1 X 10 3 M7 s, and 3.7 X 10 3 M7 s respectively. The color scale is 
logarithmic, with contour lines plotted at half order-of-magnitude intervals. 

Below some critical gas density, the Coulomb collision time 
in RIAF-type flows becomes longer than the inflow time of 
the gas, resulting in a two-temperature flow in which the ion 
plasma remains at T p ~ 10 12 K and the electron plasma cools 
to temperatures in the range T e ~ 10 10 - 10 12 K. Since elec- 
trons are more efficient radiators, the observed radiation will 
thus be determined by the properties of the electron popula- 
tion. In order to determine whether the hot plasma in the cloud 
developed a two-temperature flow, we evaluate the character- 
istic timescales for the inflow of gas fj n fl ow and Coulomb col- 
lisions fcouiomb m the vicinity of the binary. We estimate that 
for the gas cloud under consideration 



during the late inspiral, plunge, and merger of the BHs. We 
analytically estimate the characteristic timescales and emis- 
sion properties of the hot gas in Section 15.11 In Section 15.21 
we evaluate the bremsstrahlung emission from our simula- 
tions and compare it to the emission powered by accretion 
onto the two BHs (Section [5.3b . In Section IBT41 we charac- 
terize the GW signatures and highlight combined EM+GW 
signatures that may enable us to disti ngui sh binary systems 
from isolated BHs. Finally, in Section lBTBl we discuss the ob- 
servability of coalescences in the gas cloud scenario. 

5.1. Characteristic timescales and properties of the gas 

We compare the physical properties of the hot plasma in 
our simulations with the hot flows that are expected to arise in 
a fraction of AGNs as radiatively inefficient accretion flows 
(RIAFs). From these characteristic timescales, we find that 
the gas flow in our simulations in some respects resembles 
RIAFs. We do not, however, carry out a self-consistent sim- 
ulation of the multi-component plasma or radiative transfer 
in such flows as such a detailed treatment would be beyond 
the scope of this study. The main property of radiatively in- 
efficient flows is that very little energy generated by accretion 
and turbulent stresses is radiated away, being instead stored 
as th ermal energy in the gas; this gives rise to a very ho t 
flow (llchimarul 19771: iRees et al.ll982LlNarayan and Yill994l) . 
Since thermal pressure forces within the gas are significant, 
the hot accretion flows are expected to be geometrically thick. 
Given high thermal velocities, the inflow speeds in the flow 
are comparable to the speed of sound and the orbital velocity 
that a test particle would have at a given radius. 



finflow~0.4hr 



R 

Tom 



0.3c 



^Coulomb 



n ere, 



;0.5hr 



10-' 'gem- 3 



0.3c 



10'°K 



(15) 



(16) 



The size of the region under consideration here, R, is com- 
parable to the Bondi radius of gravitational influence of the 
BH binary, n rj n p rj n e w p/m p is the number density of the 
gas, c s is the speed of sound, p is the characteristic density 
of the cloud, and a w 0.3 Z 2 e 4 /{kT e ) 2 is the cross-section for 
Coulomb scattering of an electron with kinetic energy ~ kT e 
on a more massive ion, evaluated for a population of electrons 
with 7; ~ 10'° K. 

Because fcouiomb ^ Wow for the gas considered here, we 
conclude that the gas cloud region can be described as a two 
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FIG. 3. — Same as in Fig.ffjbut for the run G3. From upper left to lower right the snapshots correspond to ; = 3.1 X 10 3 M 7 s,6.2x 10 3 M 7 s,9.2x 10 3 M7S, 
1.2 X 10 4 M 7 s, 1.5 X 10 4 M 7 s, and 1.8 X 10 4 M 7 s respectively. 



temperature flow with T e < T p . However, given that the two 
timescales are not far apart in the center of our cloud, the elec- 
trons and ion plasma may remain weakly coupled and, via 
occasional scattering with ions, electrons near the BHs may 
largely preserve the thermal energy distribution. Thus, we 
estimate the bremss trahlung luminosity from a th ermal distri- 
bution of electrons dRvbicki and Lig htman 1 986b 



smaller 6 : 

^synchro ' 



x 10 36 ergs _1 
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10M 
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(18) 
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ibrem~4x 10 44 ergs 1 



p 



Mn 



10 10 K 



1/2 



10- n gcm- 



l+4.4x 



10M 



10 10 K 



(17) 



5.4 



where subscript "5.4" indicates that a numerical factor in the 
square brackets is normalized to 5.4. For comparison, the Ed- 
dington luminosity of the system is Lsdd ~ 1-3 x 10 45 ergs _1 . 
We therefore estimate that thermal bremsstrahlung is the dom- 
inant emission mechanism of the hot flow, as the synchrotron 
radiation and inverse Compton scattering are comparably 



where the relativistic factor f3 = v/c and Lorentz factor W have 
been evaluated for v/c sa 0.3. Here L so f t is a supply of low 
energy photons transported from the edge of the RIAF, a dis- 
tance of R tmn away. 

Note that the luminosity of the synchrotron emission re- 
mains below that of the bremsstrahlung radiation unless the 
magnetic field strength is close to the equipartition value, 
which in our case is B e quip ~ 10 5 G. Because our simulations 
are purely hydrodynamical, we have no means of constrain- 
ing the magnetic field strength and its dynamics around the 
coalescing binary. Thus in estimating L S ynchro> we nave as- 
sumed that magnetic fields are sufficiently weak and that the 
synchrotron component is not the dominant one. 

Similarly, in order for the inverse Compton luminosity to be 
significant, a supply of soft, lower energy photons, presum- 

6 Synchrotron and inverse Com pton scattering luminosities w ere derived 
from the standard expressions; see IRvbicki andL ightman 1983), for exam- 
ple. 
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FIG. 4. — Snapshots of the Mach number, /3., = v/c s , measured in the orbital plane of the binary in run G2, shown at the same times as in Fig. [2] Only values 
/3 S > 1 are shown. The color scale is linear ranging from light (blue) for ft = 1 to dark (red) for /3 = 5, with the black contour line marking the location of S = 1. 



ably produced externally to the hot binary accretion flow, is 
required. In estimating Lie, we assumed that this soft photon 
component is produced at large radii, where radiative cooling 
is efficient and the geometrically thick flow described here 
transitions into a geometrically thin accretion disk or a lower 
temperature ambient medium. In the case of Sgr A* for ex- 
ample, observations indicate that a radiatively inefficient ac- 
cretion flow extends into the ambient medium out to distances 
fit™ ~ 10 5 M away from the center ( IOuataertll2003l) . The es- 
timate for Lie obtained with this value of the transition ra- 
dius implies that inverse Compton scattering is a very ineffi- 
cient process even if a generous supply of low energy photons 
is available from the ambient medium, parametrized here in 
terms of the luminosity L so ft. 

Note that a s ubsequent paper to the current work 
dFarris et al.ll2.Q09l) suggests that the spinning BHs could am- 
plify the magnetic fields in their vicinity by several orders of 
magnitude to nearly the equipartition value. In the context of 
a scenario considered bv lFarris et alj d2009l) . where the mag- 
netic field strength is about an order of magnitude less than 
at equipartition (i.e., B ~ 10 4 G given the properties of our 
gas cloud), the estimate for luminosity of synchrotron radia- 
tion becomes L S ynchro ~ 10 45 ergs _1 (see equation [TFli. In this 
scenario, not only is L sy nchi- comparable to Lb rem , but the syn- 
chrotron radiation could feed the inverse Compton luminosity 
of comparable magnitude (equation \\% by providing an im- 
mediate source of soft photons. In the present work, however, 
we have no means of judging the strength of the magnetic 
field in the immediate vicinity of the binary and we focus 
most of our discussion on bremsstrahlung emission, keeping 
in mind the potential importance of the other two emission 



mechanisms. 

The corresponding cooling timescale of the plasma at T p ~ 
10 12 K due to bremsstrahlung radiation from T e ~ 10 10 K elec- 
trons is 



fcooi ~ 8 hr 



lO-Ugcm" 3 



10 12 K/ V 10 10 K 



-1/2 



(20) 



Note that f cool > fcouiomb > ^inflow implies that the hot gas 
plunges into the BHs before it had a chance to radiatively cool 
and settle into an accretion disk, thus justifying the initial as- 
sumption of a hot, geometrically thick gas cloud and implying 
that radiative cooling of the cloud can indeed be neglected in 
this case. It is also worth noting that the expressions presented 
in this section scale with density and thus can be used to es- 
timate luminosity components from a lower density plasma 
than that considered here. The scaling relations however 
break down for plasma densities higher than 10 _I1 gcm" 3 be- 
cause in this regime Coulomb collisions become sufficiently 
frequent as to thermalize the electrons and produce a single 
temperature plasma flow. Once the plasma flow of electrons 
and ions is fully thermally coupled, it can cool efficiently via 
electron-emitted radiation, yielding an evolution more similar 
to the accretion disk scenario. 



5.2. Bremsstrahlung emission from the gas 

We now evaluate the characteristic luminosity arising from 
the gas near the BBH in our simulations. At a given time 
step during the simulations, we calculate the local thermal 
bremsstrahlung emissivity £brem (i-e., luminosity per unit vol- 
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FIG. 5. — Same as in Fig.f4]but for the run G3. 



To obtain the bremsstrahlung luminosity Lbi-em, we inte- 
grate £brem over a spherical volume of radius 10M, enclos- 
ing the binary orbit but excluding the volumes inside the 
AHs. By repeating this procedure at regular time intervals 
during the simulations, we construct a light curve for thermal 
bremsstrahlung. Figure [8] shows the light curves for all four 
runs. In this and subsequent figures, t = marks the time when 
we first find the AH of the final BH. We estimate that the ac- 
tual merger, i.e. first appearance of a common AH, takes place 
within 5 x 10 2 Mj s before this. 

We first discuss the light curves in Figure [8] represented 
by solid lines. During the early part of the inspiral, t - 
5 x 10 3 M7S, the luminosities in all systems decay from the 
Eddington luminosity level to few x 10 44 ergs -1 . At about 
f ~5x 10 3 M 7 s before the merger, when the binary enters 
the final plunge, the luminosity starts increasing, leading to a 
broad flare lasting until the binary merges. The flare reaches 
its maximum during -5 x 10 2 Mj s t 0s, a period where 
we estimate the AHs of the two BHs merge. We measure the 
maximum luminosity of the flare in GO, Gl, G2 and G3 run to 
be Lbrem = { 15.5, 7.56,4.96, 13.7} x lO^ergs" 1 , respectively, 
and estimate the error associated with these measurements to 
be between 5 and 10 %. Our measurements give a lower limit 
with errors arising because the volume integrated luminosity 
depends on the formation and geometry of the "common," fi- 
nal AH, whose early distorted (peanut-like) shape is difficult 
to localize immediately after its formation. As mentioned be- 
fore, we estimate that the common, final AH forms during 
-5 x 10 2 M7 s ?it £ 0s. For this reason, we do not include in 
Fig.[8]luminosities during this time interval. 

Another characteristic feature in the luminosity curves de- 



ume) from lRvbicki and Lightmanl dl986l) as 

/ \ 2 

£brem = 2.8 x 10 4 ergs~ cm 



10 1() K 



1/2 



l+4.4x 



lO-Ugcm" 3 
10 10 K 



, (21) 

5.4 



where we assumed that T e = 10 2 T p , with T p computed from 
the internal energy e in our simulations as 

*^L=(T-l)e. (22) 
m p c L 

Figures|6]and|7]show snapshots of the bremsstrahlung emis- 
sivity £brem f° r the gas located in the plane of a binary in runs 
G2 and G3. The time sequences of the snapshots correspond 
to those in Figs. [2] -[5] for the density and Mach number. Not 
surprisingly, the highest bremsstrahlung emissivity is found in 
the regions of high density, namely at the wakes trailing be- 
hind the BHs and, at later times, around the central bar of gas 
between the BHs. By comparing Figure [6] (emissivity) with 
Figure [4] (Mach number) in runs G2 and G3, it is evident that 
regions of high emissivity are closely associated with regions 
where shocks arise and expand into the surrounding cloud. 
In the late stages of the merger, most of the emission is pro- 
duced by the bridge of shocked gas between the BHs. This 
region decreases in size as the BHs approach coalescence and 
is rapidly swallowed by the BHs when they merge to form a 
single, final BH. As we will discuss next, this gives rise to a 
characteristic variability in the bremsstrahlung light emitted 
by the gas. 
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FIG. 6. — Snapshots of the bremsstrahlung emissivity, in units of [ergs cm ], shown in the orbital plane of the binary for run G2. The time sequence of 
snapshots corresponds to Figuresf2]and|4] The color scheme is logarithmic and contours are plotted at half order-of-magnitude intervals. 



picted in Figure [8] is a sudden drop that occurs soon after the 
BHs have merged. This feature is most dramatic in run Gl 
where the luminosity decreases by nearly 3 orders of magni- 
tude. A sudden decrease can be attributed to the disappear- 
ance of the dynamic region of high emissivity between the 
two BHs, which is rapidly swallowed by the BHs in the pro- 
cess of coalescence. 

In all cases, the luminosity eventually decays exponentially. 
The luminosity in runs GO and Gl decay in this manner imme- 
diately after merger, but runs G2 and G3 exhibit an additional 
variability before the exponential decay. This exponential de- 
cay is due to the steady state accretion of the left-over gas 
surrounding the final BH. 

Since Lbrem oc p 2 , any density variations within the emitting 
volume are mirrored in the calculated luminosity curve. The 
variability following the luminosity drop in runs G2 and G3 
is especially interesting since it results in a more gradual drop 
and also highlights some differences in the evolution of the 
gas in these two runs with respect to runs GO and Gl. Specif- 
ically, in run G2, which is the longest of the four runs, the 
spiral wakes trailing behind the BHs have sufficient time to 
interact with each other, and, as a result, the final BH finds 
itself embedded in a highly turbulent medium. Because tur- 
bulence heats the gas and creates density inhomogeneities in 
the cloud, this gives rise to the variability observed in the af- 
termath of coalescence. Similarly, in run G3, the asymmetry 
of the system seeded by the prograde and retrograde spinning 
BHs, leads to interactions of the spiral wakes early in the sim- 
ulation and gives rise to turbulence and excess luminosity af- 
ter the coalescence. Since the medium surrounding the binary 
can be expected to be turbulent in realistic cases, we argue that 



realistic light curves may resemble runs G2 and G3 more than 
cases GO and Gl. Namely, we expect that turbulence would 
arise naturally in simulated systems that have been evolved 
for a longer period of time before the coalescence. 

We now highlight additional features of bremsstrahlung 
light curves that arise when relativistic beaming and Doppler 
boosting are included in the calculation of luminosity (see the 
curves plotted with dashed lines in Fig. [8j- F° r simplicity, we 
neglect relativistic bending of photon trajectories and gravita- 
tional redshift of photons in the potential well of the binary. 
We include the special relativistic Doppler effect by multiply- 
ing the broadband bremsstrahlung emissivity £brem with the 
factor D A = (W(l - /3 cos(f?)))~ 4 where 9 is the angle between 
the line-of-sight to the observer and the velocity vector of the 
gas. It follows that, depending on the position of an observer 
relative to the orientation of the binary, the changing beaming 
pattern of the orbiting binary surrounded by emitting gas can 
potentially give rise to modulations in the observed luminos- 
ity of the system. To judge the importance of this effect, we 
evaluate the bremsstrahlung luminosity for a configuration in 
which the modulations are expected to be largest by placing a 
fiducial observer in the plane of the binary at infinity. In this 
way, the observer is placed directly in the path of the sweeping 
emission beams associated with the two BHs and can sample 
both the minimum and maximum luminosity of the system, 
depending on its orbital phase. As shown in Fig.[8j see dashed 
lines, quasi-periodic oscillations in luminosity indeed arise in 
runs GO, Gl, and G2 prior to coalescence. The amplitude of 
the variations in luminosity is approximately a factor of 2 be- 
tween subsequent peaks and troughs, oscillating about a value 



Binary Black Hole Mergers and their Electromagnetic Signature 



11 




-10 S 5 10 -10 -S ^ 10 -10 

x 1,5 W 12 an X 1,5 10 "12 am 

FIG. 7. — Same as in Fig.[6]but for the run G3. The time sequence of snapshots corresponds to Figures [3]and[5] 



-3 5- 
x l.S 10" 12 on 



lower than the unboosted luminosity as boosting into the gas 
frame reduces the energy of the emitted photon by a factor 
W~ l . Furthermore, it is the beamed emission from the shocks, 
launched by the BHs into the surrounding layers of the cloud, 
that contribute the most to the modulation of luminosity. This 
argument is strengthened when oscillations in runs GO, Gl 
and G2 are compared with those in the asymmetric case, G3, 
where the binary does not form a stable set of density wakes. 
As a result, the luminosity oscillations in run G3 are much 
weaker (not easily discernible in Fig. [8]) and appear at roughly 
half the frequency seen in the other cases. 

The discussed variability stands a chance of being seen by a 
distant observer, as long as photons emitted in the vicinity of 
the BHs are not absorbed within the cloud and the surround- 
ing medium or reprocessed in such a way that the variability 
signature is lost. To estimate these effects, we consider sepa- 
rately the optical depth in the portion of the gas which serves 
as the source of bremsstrahlung photons and the remainder 
of the gas cloud. Assuming that the electron population at a 
temperature T e ~ 10 10 K will give rise to a range of photon 
energies up to kT e ~ 1 MeV, we estimate the optical depth for 
Compton scattering as rcompton ~ najR. For simplicity we 
replaced the cross-section for Compton scattering by that for 
Thomson scattering and neglect a factor of few discrepancy 
between the two that arises for the highest energy photons. 
Within the emitting portion of the cloud, approximately the 
central 10 M or 10~ 6 pc of the system, the number density of 
the gas cloud ranges from 10 12 - 10 15 cm -3 resulting in an op- 
tical depth of the order 10 - 10 4 . By contrast, the number den- 
sities outside this region drop to order 10 n cm" 3 and below, 
implying an optical depth of order unity and smaller. There- 



fore, if the photon escapes the region of emission it is likely 
that the variability in luminosity can be seen by an observer. 

5.3. Accretion onto the black holes 

In this section we discuss the properties and structure of the 
accretion flows modeled in our simulations. Fig. [9] shows the 
accretion rates measured across the apparent horizons of the 
BHs. The dash-dot lines at t < s mark the accretion rates 
onto the individual BHs before the merger. Since the accre- 
tion rates of the inspiraling BHs in runs GO, Gl and G2 are 
identical, only one curve is visible. In run G3, the dash-dot 
line at t < s marks the accretion rate for the prograde spin- 
ning BH and the dotted line that of the retrograde spinning BH 
(prograde direction is defined as a spin parallel to the orbital 
angular momentum). In all runs, the solid line represents the 
sum of the accretion rates of the two BHs before the merger 
and the accretion rate onto the final BH after the merger. 

It is prudent here to separate any effects on the accretion 
rate that stem from our simple assumptions about the initial 
structure of the gas cloud from physical effects that arise as 
a consequence of the binary dynamics. To do so we consider 
the evolution of the same gas cloud, replacing the BBH by a 
single stationary BH with a mass of IM ~ 10 7 M Q and spin 
a/M = 0.62. This choice of mass and spin approximates the 
final BH in the GO case. We overlay the accretion rate for this 
single BH case and the binary runs in Fig.[l0]leaving the time 
axes unshifted such that the starting moment of each simula- 
tion coincides (not the moment of coalescence as in previous 
figures), again omitting the initial transients. In all cases, the 
average gas density in the center of the gas cloud decreases 
exponentially as the gas cloud of finite size is swallowed by 



12 



Bode, Haas, Bogdanovic, Laguna & Shoemaker 




M [Mq yr" 1 ] 



FIG. 8. — Thermal bremsstrahlung luminosity as a function of time without 
(solid line) and with relativistic beaming (dashed line), calculated within a 
sphere of radius 10M about the center of mass of the binary but excluding 
the region within the AHs. Effect of beaming is calculated for an observer 
placed at infinity, in the orbital plane of the binary. Merger occurs at t ~ 0. 



the BHs. This decline is a consequence of our choice of ini- 
tial conditions and should not be regarded as a prediction of 
the simulations, as in reality gas may be continuously sup- 
plied to the BH and a leveled accretion rate maintained over 
longer periods of time. On the other hand, the excess variabil- 
ity in accretion rate noticeable in the four binary scenarios is 
a consequence of the binary motion which stirs the gas, caus- 
ing turbulence. It is this variability that we regard as a true 
signature of the binary. 

The two orbiting BHs accrete from a hot, turbulent flow in 
a Bondi-like fashion. We measure the total accretion rate of 
the BHs in the simulations to be in the range O.2-2M0yr -1 , 
in good agreement with the analytic expectation for Bondi ac- 
cretion rate of this system 



M B «O.84M yr" 



0.3c 



lO-Ugcm" 3 



Mi 



(23) 



where we assumed that the relative velocity between the gas 
and the BH is comparable to c s . At this rate of accretion, the 
gas residing in the nuclear region of size < 0.01 pc and mass 
1 % of the BH's mass will be accreted in ~ 10 5 yr. In order for 
the gas flow to persist uninterrupted on longer time scales, gas 
needs to be supplied from larger radii, either via an accretion 
disk or by accretion from the interstellar medium. 

A time scale of ~ 10 5 yr is, in principle, long enough for 
the angular momentum of the accreted gas to have an ef- 
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FIG. 9. — Accretion rates across the AHs. Dash-dot lines mark the accretion 
rates onto individual BHs before the merger, t < s. Accretion curves of pre- 
merger BHs in runs GO, Gl, and G2 overlap due to identical accretion rates. 
In G3, the dash-dot line marks the accretion rate onto the prograde-spinning 
BH and the dotted line marks the accretion onto the retrograde-spinning BH. 
In all cases solid line shows the total accretion rate of the system onto the 
AHs. Merger occurs at t ~ 0. 
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FIG. 10. — Accretion rate of gas onto the single spinning BH (solid line) 
and the total (sum of individual BH) accretion rates for the binary runs G0-G3 
(dashed lines), again omitting the initial transients. 



feet on the orientation of the BH spins 7 . Such coupling be- 

7 Note that in the reference frame co-rotating with the binary the gas has a 
rotational velocity component, even though its velocity distribution is initially 
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tween the angular momentum of the accreted gas and the 
BH rotation can lead to precession of their spin axes and 
possibly a partial alignment with the orbital angular momen- 
tum of the binary. This effect, however, is expected to be 
more efficient in accretion disks of small and moderate thick- 



ness (Nataraian and Prina 


e 1998; Nataraian and Armitage 


19991 iKineetal. 


2005; 


Lodato and Prinde 


120061 12007b 


Fragile et al. 2007 


; iBoadanovic et al.l 12007 


; iPereso et alj 


20091 iDotti et al.l 20091). Since the time scale over which we 



follow the binary and gas in our simulations is much shorter, 
we capture no spin alignment effects in our study and the BBH 
dynamics in all of our simulations is essentially indistinguish- 
able from the equivalent vacuum case. 
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FIG. 1 1 . — Velocity field of the gas plotted in the vicinity of the prograde- 
BH (left) and retrograde-spinning BH (right) ~ 3.3 X 10 3 M7 s before merger, 
a point at which the accretion rates differ by a factor of ~ 2. 

Since the BHs in our simulations accrete in a Bondi-like 
fashion, the expectation is that they draw matter from a re- 
gion enclosed within their Bondi radius, Rg = GM/c 2 s , which 
approximately delineates the zone of gravitational influence 
of a BH. While the Bondi radius does not explicitly depend 
on the spin of a BH, we find from run G3 that the spin ori- 
entation of the merging BHs does affect their accretion rates. 
In G3, the rates of accretion of the two BHs start diverging 
from each other at approximately ~ 6.4 x 10 3 M7 s before the 
coalescence and reach a factor of ~ 2 difference at the plunge. 
This difference arises as a combination of the frame dragging 
and geometry of the density wakes in the vicinity of the BHs. 
The panels in Fig. QT| highlight the differences in the velocity 
field of the gas around the prograde- and retrograde-spinning 
BHs at a time, ~3.3x 10 3 M7 s before merger, when the dif- 
ference in accretion rates is highest. Due to the orbital motion 
of the binary, the relative motion of the gas flow in Fig. [TTI 
appears to be from right to left at the position of the prograde- 
spinning BH and from left to right for the retrograde-spinning 
BH. A fraction of the gas near each BH tends to co-rotate 
with the spin of the BH. For both BHs the effect of frame 
dragging is such that the relative velocity of the gas with re- 
spect to the BH is effectively lower just below the BH and 
higher just above it, when compared to the mean velocity of 
the flow. This low velocity spot is, according to Eq.|23] where 
most of the accretion is favored to occur across the AHs of 
the BHs. In the case of the prograde-spinning BH, some of 
the gas plunges directly into the BH and some rotates around 
and is fed to it through the trailing wake. In the case of the 
retrograde-spinning BH, the gas is compressed into the front 
of the density wake which is then less effective at channeling 
the gas into the BH, leading to the lower accretion rate visi- 
ble in Fig. [9] Nevertheless, given the relatively weak implicit 

chaotic in the frame of the observer. 
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FIG. 12. — Maximum Lorentz factor of the gas outside the AHs. Post- 
merger values correspond to a steady Bondi-like accretion. Merger occurs at 
t ~0. 

dependence of the Bondi accretion rate on the BH spin, it is 
unlikely that one would be able to infer the spin magnitude or 
orientation in this mode of accretion based on the accretion 
powered luminosity curve alone. 

The discontinuities in the accretion rates observed in Fig. [9] 
at the time of merger occur due to the error in locating the 
common, initially highly deformed, AH of the final BH. This 
is the same source of error that we discussed in the context of 
the maximum bremsstrahlung luminosity measured at merger 
time. Notice that for t > 0, the exponential decay in the 
bremsstrahlung luminosity is also mirrored here as an expo- 
nential decay in the accretion rate of the final BH. Similarly, 
the post-merger variability present in runs G2 and G3 (Fig. [8} 
is also repeated in the accretion rates. As in the case of the 
luminosity, this behavior can be explained by the more tur- 
bulent flows of runs G2 and G3. It is also evident that the 
accretion rates measured in runs Gl and G2 fall below that in 
run GO. Since the orbital hang-ups of runs Gl and G2 result 
in a longer inspiral, the two BHs have more time to deplete 
the surrounding cloud of gas, leading to a lower density of 
gas in the vicinity of the binary and consequently lower ac- 
cretion rates at later times. This is also consistent with the 
behavior of the bremsstrahlung luminosity curves which are, 
during the post-merger phase of exponential decay, lower for 
the spinning-BH cases which exhibit a hang-up, Gl and G2, 
than for the non-spinning case GO. While details of the expo- 
nential decay phase can be dependent on our choice of initial 
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conditions and the structure of the cloud, the result that sys- 
tems with a higher net spin magnitude may appear dimmer 
with respect to their non-spinning counterparts is likely to be 
real. Without a knowledge of the gas density in the vicinity 
of SMBH binaries, though, it would be hard to infer the net 
value of the spin for a given system based solely on the in- 
ferred accretion curve. 

We now briefly discuss our results in the context of ear- 
lier works that con sidered binary coalescen ce in non-vacuum, 
specifically that of Ivan Meter et al.l (12009b . We compare the 
global maximum Lorentz factor of the gas calculated from 
our sim ulations (Fig.lT2~b with that in Fig. 2 of I van Meter et alJ 
(120091) . corresponding to their case where test particles have a 
random and isotropic velocity distribution. 8 Since it is likely 
that the maximum Lorentz factor is found in the immediate 
vicinity of the BH horizons, in Fig. [12] we do not include the 
maximum Lorentz factor during -5 x 10 2 M7S t 0s be- 
cause of, as mentioned before, the uncertainty in determining 
the region outside the final AH. 

We find that the maximum Lorentz factors measured for 
the non-spinning and spinning binary cases during the inspi- 
ral phase appear in good qualitative agreement in our respec- 
tive works. In our simulations, the ramp-up in the Lorentz 
factor during the inspiral and up to coalescence is reflected in 
the broad peak in bremsstrahlung luminosity for all runs, thus 
confirming that the increase in bremsstrahlung luminosity is 
due to shocks that arise in the vicinity of the binary. More- 
over, we confirm that the spinning systems (Gl and G2 cases), 
on average, have a higher maximum Lorentz factor than the 
non-rotating GO case. This behavior can be explained as fol- 
lows: in the case of spinning BHs, the fluid can travel deeper 
into the potential well of a BH without being accreted; as a 
consequence, the fluid can emerge with a higher kinetic en- 
ergy compared to the non-spinning case. This interpretation 
is correct as long as the fluid moves freely along a geodesic 
and before it encounters pressure forces from the surround- 
ing medium. We, however, do not find that overall higher 
Lorentz factors in the merger of rotating BHs result in higher 
luminosities. This is because the spinning BHs tend to de- 
plete their surroundings of gas more readily during the orbital 
hang-up phase and, in suc h a way , suppress the luminosity. 

Unlike Ivan Meter et al.l (120091) . we do not see a stronger 
pre-merger spike in the maximum Lorentz factor for the spin- 
ning BH cases. We think that this difference can be explained 
in the context of gas pressur e forces. Namely , in the test parti- 
cle treatment considered by Ivan Meter et al.l d2009l) . pressure 
gradients are not taken into account and each particle moves 
along a geodesic with an unaltered velocity. In our work, pres- 
sure forces from the surrounding gas act to modify the veloc- 
ity of the fluid and, in the specific case described above, result 
in a lower Lorentz factor. This can be seen in Figs. Qj] and 
[T4l which illustrate the dynamics of the gas in the immediate 
vicinity of the binary. During the inspiral, the motion of the 
binary leads to an ejection of gas which collides with the gas 
falling inwards, towards the binary. Because the ejected gas is 
unable to proceed further on its outward radial trajectory and 
is still within the gravitational influence of the binary, it even- 
tually falls back towards the BHs. In run G2, the gas is being 
ejected out by the binary on several occasions, giving rise to 
a turbulent motion in the gas cloud. In run G3 the dynamics 

8 Note that we evaluate the Loren tz factor given the velocity of a parcel 
of fluid, while van Meter et al. (2009) evaluate the collisional Lorentz factor 
from relative velocities of converging particles. 



of the flow seeded by the prograde- and retrograde-spinning 
BHs results in a lopsided ejection-infall event. In both cases 
this dynamically steered turbulence persists after the merger, 
leading to the post-merger variability in the bremsstrahlung 
luminosity and accretion rates. 

It is also worth noting that we find somewhat higher maxi- 
mum Lorentz factors in th e post-merger phase co mpared to 
the equivalent scenario in Ivan Meter et al.l d2009l) . bee ause 
our Lorentz factors reflect the bulk inflow of the gas into the 
merged BH, after the flow settles into a state of steady Bondi- 
type accretion. 

5.4. Gravitational Waves 

We now discuss properties of the GW signatures associated 
with our simulated systems as well as correlated EM+GW 
signatures. The gravitational waveforms from our simula- 
tions are almost identical to those of vacuum simulations, with 
point-wise differences of < 10~ 2 % for all the cases we con- 
sidered. This is because the mass of the gas that is gravita- 
tionally bound to our binaries is many orders of magnitude 
smaller than the mass of the BBH. 

As mentioned before, detection of correlated EM+GW os- 
cillations from the same object would be a smoking gun for 
a BBH system on the way to coalescence and would directly 
link a detected GW source to its EM counterpart. The os- 
cillations observed from the bremsstrahlung's relativistically 
beamed light are directly tied to the orbital dynamics of the 
binary and thus are also correlated with the frequency of the 
GWs. In order to compare the characteristic variability of 
the relat ivistically beamed bremsstrahlung light curve (Sec- 
tion l5.2b with that of the GWs, we overlay in Fig.[l5]the light 
curves of the beamed bremsstrahlung emission and the mag- 
nitude of the gravitational waveforms for the two prograde- 
spinning BH cases (Gl and G2) during the inspiral phase. 
Vertical lines have been drawn to mark the peaks of the real 
part of the waveform. Clearly the frequency of oscillations in 
the light curve closely matches the GW frequency. No sim- 
ilar variability is visibly present in the bremsstrahlung light 
curve in run G3 since the asymmetry in spins inhibits the for- 
mation of the strong, symmetric density wakes. As most BBH 
systems in the universe are expected to be unequal-mass bina- 
ries, and thus have some inherent asymmetry, the likelihood 
of observing such correlated oscillations therefore depends on 
the extent to which asymmetries can modify the variability as 
well as the dominant mechanism that powers the emission. 

Another form of EM variability associated with BH co- 
alescences has been predicted to arise due to the effect of 
GW mass loss. If present, this variability is expected to 
be most prominent in geometrically thin circumbinary accre- 
tion disks in which transient shocks "light up" the gas disk 
in respo nse to a changed gravitat ional potential of the cen- 
tral BH (iBode and Phinneyll2007l) . It has been shown that 
transient shocks are absent in moderately thick, hotter disks, 
where H / R > / (/ is the fractional mass loss due to emission 
of GW). In such flows the velocity of radial motions that arise 
in the gas due to perturbations in the gravitational potential is 
low with respect t o the speed of sound and thus incapable of 
producing shocks (lO'Neill et alj|2009h . With fractional mass 
losses at the level of 3.6%, 5.0%, 6.3%, and 3.6% in runs G0- 
G3, respectively, we find no subsequent EM variability in our 
simulated systems due to this effect. This is consistent with 
the expectations based on studies described above since our 
gas cloud scenario also falls in the class of hot and geometri- 
cally thick accretion flows. 
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FIG. 13. — Snapshots of the Lorentz factor measured in the plane of the binary in run G2. Outgoing flows are shown in red and ingoing flows in blue (color 
online only), with the black contour line separating ingoing and outgoing flows. The gas between the BHs is flowing radially outwards. The snapshots shown are 
coincident with those depicted in Figures[2]S and|6] 



Shocks and EM variability associated with the gas have 
also been predicte d to arise from the GW recoil of the post- 
coalescence BH dLippai et al.l 120081: IShields and Bonning 
20081: iMegevand etaj . 2009; Corrales et al.ll2009HRossi et al. 
2009tlSchnittman and Krolikl2008l) . From our four cases only 
the remnant BH in run G3 receives a "GW kick" as a con- 
sequence of an asymmetry in the emission of GWs. In this 
case, the up-down configuration of the spin vectors of the two 
BHs breaks the symmetry of the system and causes a fraction 
of linear momentum carried by GWs to be emitted preferen- 
tially in one direction. As a consequence of the conservation 
of linear momentum, the remnant BH recoils in the oppo- 
site direction with respect to the center of mass of the pre- 
coalescence binary with velocity of ~ 180kms _I . No prompt 
shocks arise in run G3 due to the GW recoil, where the evo- 
lution follows the remnant BH to a total displacement of only 
~ 0.2M from the center of mass. The absence of shocks and 
EM signatures in our simulations is again not surprising given 
the properties of the gas and the short timescales covered by 
our simulations. The imprints of the recoil on the gas could, 
however, be noticeable on much longer timescales, in which 
case a recoiling BH can produce shocks or a trail of hot, X- 
ray emi tting gas that extends out of the immediate center of a 
galaxy dDevecchi et alj|2009l) . 

5.5. Observability of coalescences 

In this section, we outline some observational strategies for 
the detection and monitoring of coalescing binaries based on 
the characteristic EM and GW signatures presented in this 
work. 



Consider a binary of mass 10 7 M Q at a redshift of z = 
1. With current expected LISA sensitivities, such a binary 
will be detectable by LISA during plunge and coalescence. 
Shortly before merger, the LISA error region could confine 
the location of suc h a binary to within a few tenths of a square 
degree on the sky dLang and Hughesll2006[) . The exact size of 
the error region depends on a number of factors, including the 
location of the object on the sky, BH masses, and spin orienta- 
tion. For instance, if the spin axes are misaligned, the orbital 
precession can further reduce the size of the error region to a 
few hundredths of a square degree. 

The LISA error box hence marks a region on the sky rel- 
evant to EM searches for counterparts. Given our assump- 
tions about the dominant emission mechanism in the vicinity 
of the BH binary it follows that Lboi ~ ^brem- The hot plasma 
cloud in our simulations is a natural source of high energy 
photons peaking in the 7-ray part of the spectrum at ener- 
gies kT e ~ 1 MeV. Because of the higher sensitivity of the 
X-ray over 7-ray surveys in terms of observed flux, in the re- 
mainder of the discussion we focus on the observational strat- 
egy in the context of X-ray observations. The X-ray emis- 
sion from our simulated sources would comprise a fraction 
of their bolometric luminosity and we estimate it by applying 
an empirical scaling in the form of a bolometric correction. 
Assuming a bolometric c orrection of Cx = 15.8, appropriate 
for low luminosity AGN (Ho 2009), we infer the X-ray lumi- 
nosity associated with our binary, Lx = L\, \/Cx ~ 10 43 erg s" 1 . 
A binary at z = 1 is at luminosity distance di = 6.6Gpc (as- 
suming a fiat universe with matter parameter il m = 0.27 and 
Hubble parameter /jq = 0.71). The observed X-ray flux from 
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this system is Fx ~ 10~ 15 erg cm -2 s -1 . According to our simu- 
lations, this level of luminosity will be observable for at least 
t > 400 M = 2 x 10 4 M 7 s in the frame of the binary or, for a 
source at z = 1, for t > 4 x 10 4 M7 s in the frame of the ob- 
server. 

Of the future X-ray observatories that may operate contem- 
porary with LISA, IXO and EXIST have planned sensitivi- 
ties sufficient to observe high luminosity obscured AGN at 
z ~ 0-2.5 and low luminosity AGN at z < 0.5. Both X-ray 
instruments have a planned field of view (FOV) close to 20' 
in diameter and a flux sensitivity limit of 10~ 15 ergcm" 2 s _1 
achievable in ~ 10 4 s of exposure time. Given the above es- 
timate for the X-ray flux, it follows that some of the brighter 
candidate AGN could be detected in as little as 1 hour of ex- 
posure. X-ray monitoring with ~ 1-hour sampling frequency 
would allow a follow-up of the quasi-periodic variability close 
to the coalescence until the point where the orbital period of 
the binary becomes < lhr/(l +z). In order not to be masked 
by a natural variability intrinsic to most AGN, the magnitude 
of quasi-periodic oscillations should exceed few tens of per- 
cent of the luminosity of an AGN. Given a LISA error region 
smaller than ~ 20', either X-ray instrument could cover it in a 
single exposure, hence allowing continuous monitoring of the 
X-ray luminosity curve. 

Note that binaries with masses lower than 10 7 M Q will 
be detected by LISA during their inspiral phase as well as 
plunge and coalescence, thus including more orbital cycles 
prior to coalescence. This may facilitate some EM counterpart 
searches since they can be triggered earlier, but in some cases 
earlier detection will imply a trade-off in strength of GW sig- 
nal. Because the GW signal is weaker during the inspiral 
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phase than during the plunge and coalescence, the LISA error 
region will be larger, and a day bef ore the coalescence it ma y 
encompass a square degree range dLangand Hughesll2006l) . 
In such a case, because the size of the LISA error region ex- 
ceeds the planned FOV of any X-ray instrument, an observa- 
tional strategy would require multiple exposure frames until 
the whole error region is tiled. In order to detect variabil- 
ity, the procedure would then need to be repeated and the 
error region scanned multiple times. Given a square degree 
error region and the 20' FOV of both future X-ray instru- 
ments, it would take about 9 exposures to cover the error re- 
gion once, before a new snapshot of the field can be acquired. 
Because lower mass binaries are also expected to be less lu- 
minous (L x ~ lO^ergs" 1 for 10 6 M Q binary), they would fall 
into a low luminosity AGN category and could only be moni- 
tored in the local universe out to a distance of ~ lOOMpc. A 
large initial error region and low luminosity make detection of 
the quasi-periodic variability a challenging prospect for lower 
mass systems. Nevertheless, a luminosity rise peaking at coa- 
lescence followed by a sudden drop-off can still be used as ro- 
bust signatures for the EM detection of these binary systems. 
Following a detection, multi-wavelength coverage of the ob- 
ject on the sky would lead to localization of the source with an 
even higher precision, at the level of arcseconds and higher, 
and would allow a study of properties of the host galaxy. 

While our discussion of observability of BH coalescences 
focuses on "LISA binaries" with masses of 10 7 M Q or less and 
the feasibility of coincident EM+GW detections, it is worth 
mentioning that a serendipitous EM detection of a BH coa- 
lescence may occur even before LISA, in the current era of 
observations. Such a detection would be most likely for more 
massive (~ 10 8 M Q ) and more luminous binary systems that 
could be "caught" by the current wide field of view observa- 
tories monitoring the high energy transient sky. Other forms 
of EM variability may also arise from mechanisms that were 
not captured by our simulations, such as radio, X-ray and 
7-ray outbursts due to reconnection and effects of magnetic 
field lines close to the binary. Specifically, if a magnetic field 
with near-equipartition strength is present in the vicinity of 
a binary the synchrotron and inverse Compton radiation may 
dominate the emission in millimeter and high energy bands, 
respectively. Whether these phenomena can give rise to char- 
acteristic signatures that would uniquely point to a BBH co- 
alescence event is a question of interest which remains to be 
studied in the context of MHD-NR calculations. 

6. CONCLUSIONS 

We presented the first fully general relativistic numerical 
hydrodynamics simulations of SMBH binaries in a gaseous 
environment through inspiral, plunge, and coalescence. The 
gaseous environment is a hot and turbulent gas cloud with 
physical properties reminiscent of accretion flows in low lu- 
minosity AGNs. The gas cloud was chosen as one of the char- 
acteristic scenarios representative of conditions in which pre- 
coalescence binaries may exist in galactic centers. Since the 
radial inflow speed in a hot gas cloud in our simulations is suf- 
ficiently high to prevent binary torques from evacuating most 
of the gas locally, accretion and interaction of the binaries 
with the gas continue uninterrupted throughout the merger. As 
a sample of the parameter space, we studied three symmetric, 
equal-mass BBH with spins of a/M = 0,0.4,0.6 aligned with 
the orbital angular momentum. In addition, we considered a 
fourth case of a binary with anti-aligned BH spins of magni- 
tude a/M = 0.4 parallel to the orbital angular momentum. The 



characteristic EM and GW signatures that arise from such in- 
teractions were the focus of this work. We summarize our 
most important results as follows: 

• Our simulations show that correlated EM-GW variabil- 
ity can occur in merging binary systems immersed in 
hot gas flows. Specifically, we found EM variability 
arising due to the effects of relativistic beaming and 
Doppler boosting modulated by the binary orbital mo- 
tion. In these systems the frequency of the EM oscil- 
lations is equal to that of the GWs and the maximum 
amplitude of variations in luminosity is a factor of ~ 2. 

• The variable EM emission in our simulations is pow- 
ered by shocks triggered by the orbiting BHs. While 
quasi-periodic variability is present in all cases con- 
sidered from inspiral through the plunge, it is not as 
pronounced in the case of a binary with asymmetric, 
prograde-retrograde spin configuration. 

• In cases where quasi-periodic variability in luminosity 
may be weak or absent, additional signatures may be 
sought for in searches for EM counterparts. Our mod- 
els indicate that, in cases where the luminosity is domi- 
nated by emission from the shocked gas, light curves 
may exhibit a gradual rise, arriving at a peak at the 
time of coalescence, followed by a sudden drop-off. 
If present, these two features are sufficiently robust 
to allow identification of an EM counterpart to a GW 
source. 

• We estimate that most massive binaries detectable in 
the LISA band may be identified in EM searches out to 
Z » 1. However, lower mass binaries and systems with 
gas densities lower than those considered here would 
fall into a class of low luminosity AGNs that could only 
be identified in the local universe. 

In summary, if coincident variability can be observed in both 
light and GWs from the same object, this signature would 
be convincing evidence for an impending BBH coalescence. 
Our results suggest some encouraging prospects for such de- 
tections. Nevertheless, given the extent of the parameter 
space involved in coalescing BBHs interacting with gas, more 
follow-up work is needed. In particular, since most of the 
SMBH binaries in the universe are expected to involve un- 
equal masses and general spin orientations, it is important 
to further explore the parameter space to investigate how 
common situations with characteristic EM signatures such as 
those in the present work may be. Furthermore, because the 
thermodynamic properties of the surrounding gas can signif- 
icantly influence the properties of EM signals, in the future 
we will also consider other scenarios for gaseous environ- 
ments around SMBHs in galactic centers, such as a rotating 
gas cloud and circumbinary accretion disk. 
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